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ABSTRACT 


A  mathematical  characterization  of  nonlinear  interpolating  spline 
curves  is  developed  through  a  variational  calculus  approach,  based  on 
the  Euler-Bernoulli  large-deflection  theory  for  the  bending  of  thin 
beams  or  elastica.  Algorithms  previously  used  for  computing  discrete 
approximations  of  nonlinear  interpolating  splines  are  discussed  and 
compared.  The  discrete  natural  cubic  interpolating  spline  is  discussed. 
An  algorithm  for  computing  discrete  natural  cubic  splines  is  given  and 
analyzed  for  discretization  error  and  computational  difficulty.  Finally, 
a  new  algorithm  together  with  its  Fortran  implementation  is  given  for 
computing  discrete  nonlinear  spline  functions. 
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1.  Introduction 


» 


For  the  past  25  year-  and  particularly  the  past  decade,  the  subject 
of  spline  functions  has  been  an  area  of  great  mathematical  interest.  The 
name  spline  comes  from  a  very  old  technique  in  drafting  in  -which  a  long 
thin  strip  of  wood,  called  a  draftsman's  spline,  is  used  to  pass  a  smooth 
curve  through  a  set  of  points  in  the  euclidean  plane.  The  points  of 
interpolation  are  called  knots  and  the  spline  is  secured  at  the  knots 
by  means  of  lead  weights  called  ducks .  The  wooden  (and  now  plastic) 
splines  and  lead  ducks  are  still  manufactured;  however,  less  expensive 
modern  drafting  tools  are  generally  used  today. 

The  mathematical  model  of  a  spline  is  a  special  case  of  the  elastic 
line,  or  elastica,  the  first  problem  of  any  importance  treated  by  the 
theory  of  elasticity.  Its  treatment  began  with  Janies  Bernoulli  in  1705, 
Daniel  Bernoulli  (1742),  Euler  (1744),  Kirchhoff,  and  many  others.  The 
history  and  theory  are  summarized  by  A.E.H.  Love  (1927);  however,  research 
papers  on  the  elastica  have  been  published  more  recently. 

Perhaps  the  simplest  way  to  characterize  a  spline  mathematically 
is  with  the  fact  that  a  spline  assumes  a  shape  which  minimizes  its 
elastic  strain  energy.  Daniel  Bernoulli  (1742)  first  suggested  that  the 
strain  energy  is  proportional  to  the  integral  of  the  square  of  the  curva¬ 
ture  taken  along  the  curve.  He  suggested  in  a  letter  to  Euler  that  the 
differential  equation  of  the  elastica  could  be  found  by  making  the  integral 
a  minimum.  Euler  (1744),  acting  on  this  suggestion,  was  able  to  find 
the  differential  equation  using  techniques  now  known  as  Calculus  of  Var¬ 
iations  and  Lagrange  multipliers,  (it  is  interesting  that  Euler  did  this 
work  when  Lagrange  was  a  small  child.) 
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Since  a  modern  development  of  Bernoulli's  discovery  could  not 
be  found  in  the  literature,  one  is  presented  here.  The  reader  is 
referred  to  Y.  C.  Fung  (1965)  for  discussions  of  the  basic  concepts 
of  solid  mechanics  used  in  the  following  discussion. 

Assume  that  the  spline  is  made  from  a  linearly  elastic  material 
with  coefficient  of  elasticity  E  .  Let  the  spline  be  initially  straight 
and  assume  that  every  cross  section  of  the  spline  which  is  initially 
perpendicular  to  the  major  axis  remains  plane  and  perpendicular  to  the 
principal  axis  at  all  times.  Let  ds  be  the  differential  arc  length 
along  the  neutral  axis.  When  the  beam  is  bent  into  a  curve  of  radius 
R  ,  the  length  of  a  filament,  initially  of  length  ds  and  parallel  to 
the  neutral  axis,  is  changed  by  a  factor  of  1  +  T|  /R  ,  where  1]  is  the 
distance  from  the  neutral  axis  to  the  filament  measured  in  a  direction 
away  from  the  center  of  curvature,  as  shown  in  Figure  1.1. 


Figure  1.1.  A  spline  element 


Thus,  the  strain  is  T|/R  ,  or  71  h  ,  where  h  is  the  curvature. 
If  dA  represents  the  cross-sp  :  ional  area  of  the  filament,  El^dA 
is  the  force  (tension)  acting  on  the  filament.  The  resultant  moment 


of  forces  acting  on  the  spline  at  the  cross  section  is 


M  =  J T]  *ET]k4A  ■  Ex y*Tl2dA 


By  definition,  the  moment  of  area  of  the  cross  section  is 


■/ 


1=1  T]  4A  . 


Therefore, 


M  =  EIx  . 


The  angle  through  which  the  cross  section  rotates  is  xds  .  Hence,  the 
work  required  to  bend  a  segment  ds  of  the  spline  to  a  curvature  h  is 


given  by 


Wk  =  (El/2)x2ds  . 


(1.1) 


There  is  a  longitudinal  force  acting  on  the  ends  of  the  element.  However, 
the  strain  energy  resulting  from  this  force  is  negligible.  Hence,  inte¬ 
grating  throughout  the  spline,  we  have  the  strain  energy 


-/ 


EIx2ds 


If  the  material  is  homogeneous,  and  the  spline  is  of  uniform  cross 


section,  then 


it  El  /  2 

U  =  ~~2~  /  *  ds  . 


(1.2) 


Now  let  It  denote  the  potential  energy  of  the  system.  In  the 


general  case 


/Wdv  -  ff.u.dv  -  /  t.u.c 

J  11  J  11 
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where  W  is  the  strain  energy  function,  f^  are  the  body  forces  per 
unit  volume,  t^  are  the  surface  traction  forces  per  unit  area,  and 
u..  are  the  displacements.  The  repeated  subscript  is  the  usual  tensor 
notation  for  summation  over  all  possible  values  of  the  subscript.  We 
are  assuming  that  the  f^  =  0  .  Now,  for  this  case,  the  t^  correspond 
to  the  forces  acting  upon  the  spline  at  the  nodes.  But  the  spline  is 
constrained  at  each  node,  so  that  either  =  0  ,  or  u^  is  orthogonal 
to  £.  .  Thus,  since  the  nodes  cannot  apply  torques  to  the  spline,  they 
can  do  no  work  on  it,  or 


Hence, 


V 


The  principle  of  virtual  work  states  that  when  a  body  is  in  equilibrium, 


Therefore, 


or 


W  =  0  . 


6U  =  0  , 


F 


ds  =  a  local  minimum 


0 


(1.3) 


(1.4) 


Of  all  the  continuous  curves  which  pass  in  turn  through  a  given  set 
of  ordered  points  in  the  euclidean  plane,  having  continuously  turning 
tangents,  and  piecewise  continuous  curvature  with  discontinuities  in 
curvature  permitted  on  only  a  finite  set  of  points,  those  satisfying 
(1.4)  are  referred  to  as  nonlinear  interpolating  spline  functions,  or 
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simply,  nonlinear  splines.  Birkhoff  and  de  Boor  (1965)  gave  an  example 
vhich  denies  the  existence  of  nonlinear  splines  in  general.  (See  Section 
3.2  for  a  discussion  of  this  example.)  It  is  an  open  question  whether 
uniqueness  holds  for  a  given  set  of  interpolation  points.  This  question 
is  complicated  by  the  fact  that  the  value  of  (1.4)  can  be  made  as  small 
as  desired  by  introducing  large  loops  between  supports,  which  modify  the 
topology  of  the  spline.  The  only  known  existence  results  for  nonlinear 
splines  are  given  in  Section  3.2. 

For  deflections  of  the  elastica  which  result  in  small  slopes, 
h  2  y"  and  ds  =  dx  .  With  this  well  known  linearizing  approximation. 
Equation  (1.4)  leads  to  the  linear  fourth-order  ordinary  differential 
equation  known  as  the  beam  equation.  The  corresponding  problems  of 
interpolation  have  unique  solutions  called  natural  cubic  spline  functions. 
These  solutions  satisfy  the  natural  boundary  conditions  of  zero  second 
derivatives  at  the  ends.  Existence  and  uniqueness  for  "linear"  cubic 
splines  follow  from  basic  theorems  of  linear  elasticity. 

As  one  would  intuitively  expect,  nonlinear  splines  are  invariant 
under  rigid  rotations  of  the  x-y  coordinate  system.  However,  the 
linear  cubic  splines  are  not  invariant  under  rotations  of  the  x-y  coor- 
ainate  system,  and  hence  they  are  not  well  suited  to  fitting  geometrical 
data  in  a  euclidean  plane,  or  other  data  where  rigid  rotations  of  the 
coordinate  system  make  sense.  On  the  other  hand,  nonlinear  splines  are 
not  invariant  with  respect  to  changes  in  the  y  coordinate  only,  while 
linear  splines  are  (see  Podolsky  and  Denman  (1964)). 

As  pointed  out  by  Lee  and  Forsythe  (1973),  the  term  nonlinear 
spline  has  been  used  variously  in  the  literature.  However,  in  this  case, 
the  term  indicates  that  the  spline  satisfies  nonlinear  Euler  equations 
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which  will  be  derived  in  the  next  section. 

Little  of  the  classical  work  on  the  elastica  is  directly  applicable 
to  nonlinear  spline  theory.  The  earliest  known  paper  on  nonlinear  splines 
is  that  of  Birkhoff  and  de  Boor  (1965).  Early  reports  were  written  by 
Fowler  and  Wilson  (196.?)  and  Birkhoff,  Burchard  and  Thomas  (1965). 

Methods  for  computing  nonlinear  splines  have  been  published  by  Glass  {1966), 
Larkin  (1966J,  and  Woodford  (1969).  Mehlum  (1969)}  in  his  Ph.p.  dissert¬ 
ation,  discusses  nonlinear  splines  and  presents  an  algorithm  for  approx¬ 
imating  a  nonlinear  spline  by  a  succession  of  circular  arc s.  Hosaka  (1967) 
describes  how  to  solve  (1.4)  and  generate  nonlinear  spline  functions  on 
a  differential  analyzer.  Some  of  these  algorithms  will  be  discussed  in 
Section  3- 

Methods  for  computing  nonlinear  splines  nave  applications  in  the 
design  of  aircraft,  ships  and  automobiles.  A  piece  of  sheet  ir  'tal  which 
is  bent  into  an  axially-symmetric  configuration  has  a  cross  section  which 
behaves  lile  a  nonlinear  spline. 


* 
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A  Variational  Formulation  of  the  Large-Deflection,  Small-Strain 
Problem  of  Interpolating  Eiastica 

Given  a  fixed  sequence  of  points  P  ,  P2,  . ..,  Pn,  we  define  the 
interpolating  eiastica  as  the  function  P(s),  a  function  of  arc-length 
s,  which  satisfies  the  equilibrium  conditions  for  a  thin  beam,  of  con¬ 
stant  cross  section  and  constant  linearly  elastic  properties,  passing  in  turn 
through  the  n  points  P  and  acted  on  only  by  workless  constraints 
at  the  points  P^  .  As  we  shall  see,  the  natural  boundary  conditions 
resulting  from  the  variational  formulation  will  indicate  exactly  what 
kinds  of  boundary  conditions  are  admissible  and  what  combinations  of  them 
result  in  well-posed  problems.  We  will  also  see  that  the  forces  f\  acting 
on  the  eiastica  through  the  points  p_^  must  satisfy  certain  restrictions 
which  are  entirely  consistent  with  the  assumptions  of  Section  1. 

It  is  natural  to  specify  each  P^  by  a  pair  of  Cartesian  coor¬ 
dinates  (x^y^;  however  a  more  efficient  coordinate  system  for  the 
development  of  the  problem  is  (0,s),  where  s  is  the  arc  length  and 
0  =  0(s)  is  the  angle  made  at  s  by  the  curve  P(s)  with  some  refer¬ 
ence  line.  The  curvature  of  the  eiastica  at  any  point  s  is  given  by 
h  ( s )  =  d0/  ds  . 

In  Section  1  it  was  shown  that  stable  configurations  of  the  eias¬ 
tica  correspond  to  local  minima  of  the  functional 


(2.1) 


where  l  is  the  length  of  the  eiastica  and  the  integral  is  taken  along 
the  deformed  configuration.  Note  that  l  is  allowed  to  vary  in  minimizing 
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(2.1).  It  is  equivalent  to  say 


A  /  h  ds  =  0  . 


(2.2) 


For  the  trivial  case  of  collinear  P.^  and  all  f^  =  0  ,  the  value  of 
(2.1)  for  the  equilibrium  state  is  0  .  To  find  the  function  0(s)  which 
satisfies  (2.2)  for  a  nontrivial  set  of  P  ,  the  additional  2n-2  side 


conditions 


cos9  ds  =  xi+1  -  xi  ,  i=l,  . ..,  n-1 


( 


sinq  ds  =  yi+1  -  y±  ,  i=l,  . . . ,  n-1  , 


(2.3) 


must  be  added  to  (2.2),  where  denotes  the  value  of  s  at  P^  . 

Using  Lagrange  multipliers,  the  constraints  (2.3)  can  be  added  to  (2.2) 


yielding 


n-1  ..  A+l 


r  r  2 

6  <  y  II  (h  +  l^cos  0  +  p^sin  fi)ds 
1=1  ii 

"Xi(xi+1  -  V  -  “l^i+l  '  Vj  }  =  0  • 


Now  since  certain  kinds  of  constraints  at  the  P^  allow  the  elastica  to 
slide,  the  lengths  are  allowed  to  vary  as  well  as  the  function  0  . 
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Taking  continuous  variations,  integrating  by  parts  and  rearranging  terms 
gives 


n-1  l 


*1+1 

S  /  [‘2  S  sl” s  +  cos  e] 


60ds 


i=l  li 


+  2k~6  0  -  2hT6 0, 

n  n  11 


n-1 

+  2^°*i'  - 


i=2 


-\2 


+  [(H)  +  X  cos  0"  +  ^  ,sin  0"]&£ 

v  n '  n-1  n  n-1  n  n 

-  [(h|)2  +  X  1cot  P*  +  p^sin  9+]  6^ 


+  £  [(Hi)2  *  (<)2 

i=2 


(2.4) 


+  X  .  , cos  G .  +  a.  , sin  07 

l-l  l  l-l  l 

-  X  .cos  fit  -  a. sin  0+]  6/.  =  0  , 
i  i  i  i  1 

where  the  notation  ?t  denotes  lim  ?(£.  +  *),  and  ?7  denotes  lim  S(V«) 

»0  1  c-*0 

By  the  fundamental  lemma  of  the  Calculus  of  Variations,  the  integral 
term  in  (2.4)  yields 
ri.  k 

-2  —  -  X  .sin  0  +  p.cos  6  =  0  ,  A.  <  s  <  A.  ,  (2.5) 

ds  i  l  i  i+1 

for  i  =  l,...,n-l  .  Using  (2.5),  Equation  (2.5)  can  be  integrated  for 

each  open  interval  yielding 

+  X,  M-i 

k(s)  =  -  -g-  (y-yi)  +  -g-  (x-x^  ,  i=l, . . . ,n-l  .  (2.6) 

Equation  (2.6)  also  results  from  the  static  equilibrium  moment  relation  for 
a  segment  of  spline  to  the  right  of  the  i-th  knot  (see  Figure  2.1). 
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Figure  2.1.  A  segment  of  spline  to  the  right  of  the  i-th  knot 

Hence,  the  Lagrange  multipliers  are  scaled  force  components  acting  on  the 
spline  at  . 

The  Euler  equation  can  be  deduced  by  differentiating  (2.5)  which  gives 

^  \ 

-  -f-  cos  0  -  sin  0  1  K  ,  (2.7) 

So,  if  h  4  0  , 


Hence,  for  h  /  0  , 

2 

—  —  =  -^h  +  c  ,  i.  <  s  <  ,  i  =  1, . . .  ,n-l  ,  (2.<t) 

K  ds2  x  i  1+1 

for  constants  depending,  in  general,  upon  the  interval  i  . 

The  remaining  terms  of  (2.4)  yield  the  natural  boundary  conditions. 
Since  each  of  the  Sfb  and  6^  are  arbitrary  (for  i=l,...,n),  the  knots 
of  the  spline  can  be  modeled  by  frictionlessly  rotating  sliders,  as  shown  in 
Figure  2.2. 
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The  terras 

2<bMn  -  2»X 

in  (2.4)  give  the  natural  end  conditions  of  zero  curvature.  The  terms 

n-1  -  + 

2  E  (h  -  h  )60 

i=2 

demand  that  the  curvature  be  continuous  across  each  interior  knot.  If 
the  angle  0  is  taken  with  reference  to  the  x-axis,  we  have  f  'om  (2.7) 

T(s)  =  -  ~  ,  (2.9) 

K  dsd 

where  T(s)  denotes  the  (scaled)  force  transmitted  along  the  spline. 
(Tension  is  positive.)  Thus,  the  terms 

[(h")2  *  \n  .cos  0"  +  u  .sin  0']6i 
n  n-l  n  n-.l  n  n 

-[(k*)2  +  X^cos  0;  +  i^sin  6^)6^ 

in  (2.4)  give  the  additional  natural  end  conditions  T(4*)  =  T(^n)  =  0  . 
That  is,  the  ends  of  the  spline  have  no  longitudinal  forces  applied  to 
them.  Similarly,  the  terms 

E  [(O2  -  (k+)2 
i=2 

+  Xi_1cos  0^  +  M-i_1sin  0^ 

-  X.cos  0?  -  u.sin  fit]80. 
i  l  ri  i  i 
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require1  that  the  forces  T(4^)  =  T(jtT)  (i=2 , . . . ,n-l)  ,  that  is,  the  longitud¬ 
inal  force  transmitted  by  the  spline  is  continuous  across  each  interior  knot. 


The  constants  of  integration  c^  (i=l, . . .  ,n-l)  in  (2.8)  can  now  be 
determined  for  the  natural  open  spline.  Since  T(i1)  =  k  =  0  ,  (2.8) 
evaluated  at  gives  c^  =  0  .  Both  the  curvature  and  the  longitudinal 

force  are  continuous  across  each  of  the  interior  knots,  hence  we  have 
e,.  --  0  for  i=2,...,n-l  .  Thus  we  obtain  the  Euler  equations 


TS  +  S*’  =°  *  l±<  s<  *i+1  » 


(2.10) 


for  i=l,...,n-l  ,  which,  in  view  of  (2.7),  must  be  satisfied  by  the  spline 
at  every  point  within  each  open  interval.  Equations  (2.10)  were  first 
published  by  Birkhoff  and  deBoor  (1965). 

The  above  analysis  shows  that  the  nonlinear  spline  having  friction- 
iessly  rotating  slider  constraints  corresponds  to  the  least-constraint 
interpolating  elastica.  Any  further  constraints  on  the  elastica  must  be 
consistent  with  (2.1)  and  will  produce  a  larger  value  of  (2.1)  unless  they 
are  also  consistent  with  the  nonlinear  spline.  For  example,  Hermite  con¬ 
straints  arc  consistent  with  (2.1)  since  they  imply  60^  =  0  ,  (i=l,...,n)  . 
However,  for  this  case,  curvature  is  not  necessarily  continuous  across  the 
interior  knots  and,  hence,  (2.10)  no  longer  holds.  Numerous  combinations 
of  more  restrictive  boundary  conditions  are  admissable  according  to  (2.1). 
These  include  such  things  as  fixed  ends  and  pinned  joints.  We  will  only 
be  ejn'erned  with  the  least-constraint,  natural  boundary  conditions  of  the 
nonlinear  interpolating  splines. 

Closed  nonlinear  splines  have  been  studied  by  Lee  and  Forsythe  (1^7  i). 
In  this  case,  and  y^  =  yn  .  However,  to  avoid  having  to  pres¬ 
cribe  the  total  length  of  the  spline,  the  variations  and  64^  are  taken 
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as  unequal.  Since  elements  ol'  the  curved  arc  are  inserted  into  or  deleted 


l’rom  the  spline  at  =  >  the  variations  at  t  and  must  satisfy 


6  6,  >  h,  61,  =  69  +  k  61  . 

1  11  n  n  n 


The  corresponding  natural  boundary  conditions  are 

Hl  =  Kn 
and 

(Kp2  +  2T*  =  (k“)2  +  2?~n  =  0  . 

It  follows  that  (2.10)  also  holds  for  the  case  of  closed  nonlinear  splines. 
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3-  Methods  for  Computing  Open  Nonlinear  Spline  Functions 

A  number  of  methods  for  computing  nonlinear  interpolating  splines 
have  been  published  during  the  past  six  years.  In  this  section  these 
results  are  briefly  reviewed  in  chronological  order.  Also,  two  interesting 
results  due  to  Larkin,  which  concern  the  existence  of  finite  equilibrium 
solutions,  are  treated  in  Section  3*2. 


3-1  Glass'  method 

The  first  published  method  for  computing  discrete  approximations 
to  the  nonlinear  interpolating  spline  is  that  of  Glass  (19 66). 

Glass  uses  Cartesian  coordinates  (x-y)  and  considers  the  Euler 
equation  analogous  to  (2.10): 

y(iv)  =  M/)3  +  .g9y  (3  x 

y  «i*(x  'ff 

The  knots  are  represented  by  their  coordinates  (x^y^),  i=l,2,...,n  . 
Since  the  curvature 


k  = 


[l  +  (y')2]5/2 


v 


and  ds  =  4/  1  +  (y')^  dx  ,  the  energy  given  by  (2.1)  becomes 
n-1  x,, 


E(y)  = 


y  [i+l _ tip2  . 

La  J  [1  +  (y')2]5^2 


dx  . 


(3.2) 


i=l 


Glass  first  considers  the  problem  on  one  panel  [  x^,x^],  with  the 
end  points  y^^  and  yi+1  prescribed.  He  assumes  the  end  slopes  y/ 
and  yi/+1  are  also  prescribed  and  proceeds  to  develop  a  method  to  solve 
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the  ordinary  differential  equation  (3.1),  subject  to  these  four  boundary 
conditions.  Using  a  Taylor's  series  expansion  of  (3.1)  about  an  approxi¬ 
mate  solution  y^°),  he  arrives  at  a  linear  differential  equation  in  the 
correction  to  y^  to  obtain  the  better  approximate  solution  y^  . 

Using  central  difference  operators  over  a  mesh  of  N  points  in  the  inter¬ 
val  (with  constant  mesh  size  H),  he  obtains  a  system  of  linear  difference 
equations  in  the  correction  terms  at  each  mesh  point.  This  5-band  system 
of  linear  algebraic  equations  is  solved  at  each  iteration  until  convergence. 

The  solution  of  the  fourth-order  boundary  value  problem  is  carried 
out  in  each  ponel  for  a  given  set  of  slopes  at  the  knots;  y'  ,  i=l,2,...,n  , 
written  y'.  The  functional  (3.2)  is  approximated  using  this  discrete 
solution  and  is  now  considered  to  be  E(y')  .  Glass  then  proceeds  to  a 
solution  with  an  outer  iteration  which  minimizes  E(y')  using  a  gradient 
method.  He  uses  differencing  to  obtain  partial  derivatives  and  then  rounds 
them  to  1,  0,  or  -1.  He  uses  an  unusual  heuristic  method  for  choosing 
the  step  parameter  X  for  moving  in  the  direction  -grad[E(y')]  . 

Glass  doesn't  present  the  actual  program;  however,  he  claims  that  his 
method  requires  only  about  five  gradient  searches  for  satisfactory  answers. 

He  presents  an  example  with  nine  panels  (n=10),  but  does  not  mention  the 
value  of  N  or  how  the  initial  y'  was  obtained.  For  this  example, 

Glass'  method  took  approximately  10  minutes  on  the  IBM  709^. 

In  addition  to  the  "enormous  computation  time,"  Glass  mentioned  that 
*  instabilities  occur  in  problems  where  the  spline  turns  too  rapidly. 

He  suggests  a  technique  for  overcoming  the  instabilities  at  the 
expense  of  a  substantial  increase  in  computation  time.  The  computational 
f  difficulties  increase  drastically  for  larger  numbers  of  data  points. 
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3  .2  Larkin's  method 

Larkin  (1966)  presented  perhaps  the  most  interesting  study  of  non¬ 
linear  spline  functions  in  a  rather  obscure  unpublished  report.  His  results 
and  suggested  computational  method  are  summarized  here. 

Larkin  first  shows  that  (3.1)  may  be  reduced  to 

*  =  “Js“  =  xi  V"c°s(e^l  ,  ii  <  s  <  jti+1  ,  (3.3) 

where  the  \ ^  and  r  are  constants  of  integration.  Equation  (3.3) 
is  precisely  the  energy  integral  of  the  equations  of  motion  of 
Kirchhoff's  kinetic  analogue  of  a  pendulum  given  in  Love  (1927 )>  P-  ^01, 
which  may  also  be  written  as 

k  =  Ai  cos  8  +  Bi  sin  0  ,  ii  <  s  <  Xi+1  ,  (3.*0 

where  the  and  are  the  constants  of  integration.  Equation  (3.4) 

is  easily  shown  to  satisfy  (2.10)  if  there  are  no  points  of  inflection. 

Larkin  integrates  (3.3)  to  give  x(s)  and  y(s)  in  terms  of  elliptic 
integrals.  To  obtain  equations  which  can  be  used  for  computation,  he  needs 
a  separate  treatment  for  the  case  where  the  spline  has  a  point  of  inflection. 
A  point  of  inflection  occurs  when  the  curvature  vanishes,  i.e.  where 

e  =  e.  +  rr/2  .  (3.5) 

Larkin's  complete  algorithm  is  quite  complicated  and  goes  roughly  as 
follows : 

Step  1:  Estimate  values  of  0  at  the  knots  (i.e.,  0^  ,  i=l,2, . . . ,n) . 

Step  2;  Determine  the  parameters  \ ^  and  ,  for  each  panel.  This 
requires  the  evaluation  of  three  complicated  functions,  each  involving 
elliptic  integrals,  as  well  as  several  trigonometric  functions.  Based 
upon  the  values  of  these  functions,  it  can  be  determined  whether  or  not  a 


l6 


point  of  inflection  occurs  in  the  interval,  and  if  so,  which  one  of  two 


cases  it  corresponds  to.  The  value  of  is  then  determined  by  finding 
the  zero  of  the  appropriate  complicated  transcendental  equation.  The 
value  of  \  is  then  obtained  by  the  direct  evaluation  of  a  similar 
equation. 

Step  3:  Scan  through  the  knots  to  find  the  one  at  which  the  largest  dis¬ 
continuity  in  curvature  occurs.  Replace  the  value  of  there  by  a  value 

which  makes  the  curvature  continuous.  This  requires  solving  for  the  zero 
of  another  transcendental  equation,  and  the  values  of  and  X  must 

be  recomputed  in  each  of  the  neighboring  panels.  Step  3  is  repeated  until 
the  largest  discontinuity  in  curvature  becomes  smaller  than  a  specified 
tolerance. 

The  above  algorithm  was  implemented  on  a  KDF9  computer  and  graphs  of 
a  few  solutions  are  included  in  Larkin's  report.  No  mention  is  made  of 
running  times,  or  computational  experience.  However,  it  appears  that 
the  method  is  probably  quite  slow  dve  to  the  large  number  of  transcendental 
functions  which  must  be  evaluated. 

The  most  interesting  results  in  Larkin's  report  are  contained  in 
the  following  two  theorems; 

Theorem  3 . 1  (Larkin):  A  necessary  condition  for  the  existence  of  a  finite 
equilibrium  solution  is 

|  "  0-jJ  S  ’  i=l,2 , . . .  ,n-l  .  (3>6) 

Proof ;  From  Kirchhoff's  analogy  (5.3), 

2  2 

H  =  X.  cos  ( 0  -  * . )  >  0  . 

Thus,  cos  (0  -  *j)  >  0  .  Therefore,  at  4^  and  ,  we  have 
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| ©i  -  cj  <  n/2  +  2pn  , 

and 

|6i+l  -  «i!  <"/2  +  Sqn  , 

where  p  and  q  are  integers .  The  only  cases  of  interest  are  those 
where  p  =  0  ,  and  is  evaluated  so  that  p  =  q  =  0  .  So  we  have 

I  (0i+l  "  *i^  "  ^®i  ei^i  -  I ®i+l  ’  *i I  +  l0i  €il  »  H 

or  I  0i+l  “  ®i I  —  n  *  I 

Unfortunately,  this  theorem  doesn't  provide  a  way  to  determine  a 
priori  whether  a  given  set  of  knots  has  a  finite  equilibrium  solution. 
However,  this  result  is  useful  in  an  algorithm  for  determining  when  a 
solution  is  diverging.  The  cases  where  (3.6)  is  not  satisfied  correspond 
physically  to  a  spline  which  continues  to  slide  through  the  supports  while 
reaching  lower  and  lower  energy  values.  An  example  of  this  is  given  in 
Birkhoff,  Burchard  and  Thomas  (1965).  The  knots  in  this  example  have  the 
Cartesian  coordinates  (l,0),  (2,0),  (0,2),  (0,l),  with  the  spline 
threaded  through  tuem  in  that  order.  The  spline  never  reaches  an  equil¬ 
ibrium  state,  and  continues  expanding  to  infinity.  A  lucid  explanation 
of  this  example  is  given  in  Lee  and  Forsythe  (1973). 

Theorem  3.2  (Larkin):  A  necessary  condition  for  the  existence  of  a  finite 
equilibrium  which  possesses  no  point  of  inflection  is  that  there  exists  an 
integer  m  such  that 

min  (0i}6i+1)  <  ♦i  +  mn<  max  (9i>8i+1)  ,  (3-7) 

where  is  the  inclination  of  the  straight  line  connecting  the  i-th 
and  the  (i+l)-th  knots. 
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Proof;  Roughly,  this  theorem  follows  from  the  fact  that  0  is  a  contin¬ 
uous  function  which  must  vary  raonotonically  from  0^  to  0^+1  in  the 
case  that  there  is  no  point  of  inflection.  And  (3.7)  is  untrue  only  if 
0  is  not  monotonic.  For  a  detailed  proof,  see  Larkin  (1966),  Appendix  A. 

3.3  Woodford’s  method 

The  only  published  program  for  computing  nonlinear  splines  is  an 
Algol  procedure  given  by  Woodford  (1969).  His  method  appears  to  be  con¬ 
siderably  simpler  and  faster  than  those  of  Glass  and  Larkin. 

Woodford  uses  cartesian  coordinates  and  applies  the  calculus  of 
variations  to  obtain  (3.1).  He  notes  that  (3.1)  can  be  reduced  to 

(y*)2  =  (Ajy'  +B±)(i  +  (y')2)5/2  (3.8) 

which  is  merely  another  form  of  Kirchhoffs  kinetic  analogue  given  by 

(3.3)  and  (3.4).  He  uses  (3.8)  to  obtain  an  expression  for  the  energy 
n-1 

E  =  ^  '  |Ai(yi+i  "  y*)  +  "  xi)|  »  (3*9) 

i=l 

which  is  easily  derived  using  (3.4)  and  (2.1). 

Woodford  discretizes  the  intervals  with  a  uniform  mesh  and 
starts  with  an  initial  solution  y^  provided  by  solving 


which  gives  the  usual  natural  cubic  interpolating  spline.  He  uses  the 
method  of  quasi-linearization  and  proceeds  with  an  iteration  based  on  the 
same  Taylor  series  expansion  about  the  current  solution  that  Glass  used. 
To  solve  this  linear  multi -boundary  value  problem,  Woodford  uses  a  deriv¬ 
ative  replacement  scheme  based  on  expressing  y^+1,  y^+1,  y£+1>  and  y* 
with  Taylor  series  expansions  about  the  point  y^  ,  immediately  to  the 
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.  This  results 


( iv) 

left,  and  retaining  all  terms  up  through  those  in  y 
in  a  9-band  linear  algebraic  system  of  equations  to  solve  at  each  iteration. 
The  iterations  are  terminated  when  the  energy  given  by  (3*9)  on  two  succes¬ 
sive  iterations  is  nearly  the  same. 

Woodford  gives  an  example  using  seven  data  points:  (0,0),  (1,1.9), 
(2,2.7),  (3,2.6),  (4,1.6),  (5,0.8),  (6,1.2),  and  a  mesh  size  of  0.025 
(40  points  per  panel).  Convergence  occurred  after  four  iterations.  His 
example  will  be  discussed  again  in  Section  5*  A  listing  of  a  Fortran  trans¬ 
lation  of  Woodford's  Algol  procedure  follows. 


SUBROUTINE  MEC (N,  ORD,  H,  K,  A,  B,  C,  D,  EPS,  *) 


C 

C  THIS  SUBROUTINE  COMPUTES  POINTS  ON  THE  MINIMUM-ENERGY  CURVE 
C  (SOMETIM5S  CALLED  A  NONLINEAR  SPL  WHICH  PASSES  THROUGH 
C  THE  N  DA' ’A  POINTS  WITH  ORDINATES  ORD(i),  1-1,..., N,  AND 
C  EQUIDISTANT  ABSCISSAE.  THE  SOLUTION  IS  PRODUCED  BY  A  STEP-BY- 
C  STEP  METHOD  BASED  ON  TAYLOR  SERIES.  THE  NUMBER  OF  STEPS  BETWEEN 

r  CONSECUTIVE  PAIRS  OF  DATA  POINTS  IS  K,  AND  H  IS  THE  LENGTH  OF 

C  EACH  STEP.  USUALLY  THE  NUMBER  OF  STEPS  K  SHOULD  BE  AT  LEAST  10 
C  PER  UNIT  INTERVAL.  THE  PARAMETER  EPS  CONTROLS  THE  CONVERGENCE, 

C  USUALLY  EPS  SHOULD  BE  OF  THE  ORDER  1.E-5.  THE  SOLUTION  Y  AND  ITS 
C  DERIVATIVES  ARE  STORED  IN  THE  ARRAYS  A,  B,  C,  AND  D.  A(I) 

C  IS  THE  VALUE  OF  Y  AT  THE  I-TH  POINT,  1-1 . K*(N-1)+1.  SIMILARLY 

C  B (I)  IS  THE  VALUE  OF  THE  FIRST  DERIVATIVE,  C(I)  THE  SECOND 

C  DERIVATIVE  AND  D(I)  THE  THIRD  DERIVATIVE.  SINCE  WE  ALLOW  FOR 
C  DISCONTINUITIES  IN  THE  THIRD  DERIVATIVE  AT  EACH  INTERNAL  DATA 
C  POINT,  D(I)  FOR  I  -  K+1,  2K+1 ,  ....  (N-I)K-M  IS  THE  THIRD 
C  DERIVATIVE  TO  THE  IMMEDIATE  RIGHT  OF  THE  DATA  POINT. 

C  D((N-1)*K+1)  IS  LEFT  UNDEFINED.  THE  SUBROUTINES  BANDET  AND 
C  BANDSL  MUST  BE  SUPPLIED  FOR  SOLVING  A  5-BAND  LINEAR  SYSTEM 
C  AT  EACH  ITERATION.  THE  SUBROUTINE  EXITS  WITH  A  RETURN  1  IF 
C  THE  PROCEUDRE  BANDET  FINDS  A  ZERO  PIVOT  OR  IF  THE  ALGORITHM 
C  HAS  NOT  CONVERGED  AFTER  10  ITERATIONS. 

C 

C  THIS  SUBROUTINE  IS  A  TRANSLATION  OF  AN  ALGOL  PROCEDURE  GIVEN 
C  BY  C.H.  WOODFORD  IN  "SMOOTH  CURVE  INTERPOLATION",  BIT  9  (1969), 

C  69-77. 

C 

C  MICHAEL  A.  MALCOLM 
C  COMPUTER  SCIENCE  DEPARTMENT 
C  STANFORD  UNIVERSITY 
C  AUGUST  16,  1971 
C 

REAL*4  H,  EPS,  0RD(1),  A(1),  B(1),  C(1),  D(1),  AL(1000,9), 

*  R(1000),  EN(20) ,  DEN(20) ,  AA,  BL,  CL,  DL,  X,  W,  ZZ,  XX,  XXI, 

*  XX2 ,  WW,  WWW,  D1 ,  HI,  H2,  H3,  H4,  MBAND(1000,  4),  NORM 
INTEGER  Q,  P,  INT(IOOO) 

LOGICAL  L,  Z 
C 

D1  -  1.E10 
ITO  -  0 
N1  -  N-1 
HI  «  H*K 
H2  -  . 5*H*H 
H3  =  H2*H/3 . 

H4  =  H3*H*.25 
Q  ■=  N1  *  (4*K-1 ) 

Ml  -  N1*K  +  1 
M  -  MI-1 


> 


r 


f 


♦ 


DO  10  1*1, Ml 
A(I)  *  0. 

B ( I )  -  0. 

C(I)  -  0. 

D ( I)  -  0. 

10  CONTINUE 
EN(1)  -  0. 

EN(N)  -  0. 

20  DO  30  1-1, Q 

DO  30  J-1,9 
AL(I, J)  *  0. 

30  CONTINUE 
1*5 
J  -  3 

DO  80  P*2,M 

L  =  (P/K)*K.EQ.P 
Z  *  ((P-1)/K)*K.EQ. (P-1) 
AA  -  A(P) 

BL  -  B(P) 

CL  -  C(P) 

DL  -  D(P) 

11  *  1+1 

12  =  1+2 

13  *  1+3 

14  *  5-1 

15  *  4-1 

16  -  3-1 

17  =  2-1 

IF  (Z)  GO  TO  40 
AL(I,J+I4)  =  1. 

J  =  J+1 

40  AL(I , J+I4)  -  H  +  CL*H4 

AL(I1,J+I5)  -  1.  +  CL*H3 
AL(I2,J+I6)  =  CL*H2 
IF  (L)  GO  TO  42 
AL (13,  J+I7)  *  CL*H 
42  J  *  J+1 

AL(I, J+I4)  *  H2  +  BL*H4 
AL(1 1  ,J+I5)  =  H  +  BL*H3 
AL(I2,J+I6)  *  1.  +  BL*H2 
IF  (L)  GO  TO  44 
AL(I3, J+I7)  -  H*BL 
44  J  -  J+1 

AL(I , J+I4)  -  H3  +  AA*H4 
AL(I1 , J+I5)  =  H2  +  AA*H3 
AL(I2 , J+I6)  *  H  +  AA*H2 
IF  (L)  GO  TO  46 
AL(I3, J+I7)  -  1.  +  AA*H 
46  J  =  J+1 

IF  (L)  GO  TO  48 
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48 


AL(I , J+I4)  -  -1. 

J  =  J+1 

AL (I 1 , J+I5)  -  -1. 
J  -  J+1 


IF  (P.NE.M)  AL(I2 , J+I6)  -  -1. 

J  -  J+1 

IF  (L)  GO  TO  50 
AL(I3, J+I7)  -  -1. 

50  BL  =  DL*H4 

R(I1)  -  DL*H3 

R(I2)  =  DL*H2 

IF  ( .NOT.L)  GO  TO  60 

IPORD  -  P/K  +  1 

R(I)  =  BL  +  ORD (IPORD) 

I  =  13 
J  =  J-2 
GO  TO  80 

60  IPORD  =  (P-1)/K  +  1 

R(I)  BL 

IF  (Z)  R(I)  =  R(I)  -  ORD (IPORD) 

R(I3)  =  DL*H 
I  =  1+4 
J  =  J-3 
80  CONTINUE 
AA  =  A  ( 1 ) 

CL  *  C(1) 

DL  =  D (1 ) 

AL(1 ,5)  =  H  +  CL*H4 
AL(1 ,6)  =  H3  +  AA*H4 
AL  ( 1 , 7 )  =  -1. 

AL ( 2 , 4 )  -  1.  +  CL*H3 
AL (2 , 5)  =  H2  +  AA*H3 
AL ( 2 , 7 )  =  -1. 

AL  ( 3 , 3 )  =  CL*H2 
AL ( 3 , 4 )  =  H  +  AA*H2 
AL(3, 7)  =  -1. 

AL(4,3)  =  1.  +  AA*H 
AL(4 ,2)  =  CL*H 
AL(4 , 7)  =  -1. 

R( 1 )  =  H4*DL  -  ORD ( 1 ) 

R(2)  =  H3*DL 
R(3)  =  H2*DL 
R(4)  -  H>DL 

CALL  BANDET(MBAND,  1000,  INT,  AL,  Q,  4, 
CALL  BANDSL(MBAND,  1000,  INT,  R,  AL,  Q, 
C  STORING  CURRENT  SOLUTION 
1=4 
J  =  2 

DEN(1 )  =  R(1) 

DEN(N)  =  R(Q) 


9,  &150) 
4,  9) 
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no  100  P-2,M 
X  -  R(I) 

W  -  R(I+1) 

ZZ  -  R(I+2) 

I  -  1+3 

IF  ((P/K)*K.NE.P)  I  -  1+1 

XX  -  X*X 

XXI  -  1./(1.+XX) 

XX2  -  XXI  *  XXI 

ww  -  w*w 
WWW  -  ww*w 

IF  (((P-1  j/K)*K...E.  (P-1))  GO  TO  90 
EN(J)  -  WW* (XXI **2.5) 

DEN(J)  -  X 
J  -  J+1 

90  AA  «  10.*X*W*XX1 

BL  -  ( 7 . 5*WW  +  10.*X*ZZ)*XX1  -  52 . 5*XX*WW*XX2 
CL  -  10.*W*ZZ*XX1  -  X*WWW*XX2*(40.  -  70.*XX*XX1) 

*  -  20.*XX*W*ZZ*XX2 

D(P)  -  -(2.5*WWW  +  10.*X*W*ZZ)  *  XXI  +  1 7 . 5*XX*WWW*XX2 

*  +AA*ZZ  +  BL*W  +  CL*X 
A (P)  -  AA 

B(P)  -  BL 
C(P)  -  CL 
1 00  CONTINUE 
C  TEST  FOR  CONVERGENCE 
X  -  0. 

DO  110  1-1, N1 

AA  -  (EN(I+1)  -  EN(I))/ (DEN (1+1)  -  DEN(I)) 

BL  -  EN(I)  -  AA*DEN(I) 

X  -  X+ABS (AA* (ORD (1+1 )  -  ORD(I))  +  BL*H1) 

110  CONTINUE 

ITO  -  ITO+1 
NORM  =  ABS(DI-X) 

PRINT  115,  ITO,  NORM 

115  FORMAT ( '  ITER.  NO.  ',  14,  '  NORM  -',  G10.3) 

IF  (ABS(DI-X) .LT.EPS)  GO  TO  120 
IF  (ITO .EQ. 10)  RETURN  1 
D1  -  X 
X  -  R  ( 1 ) 

B(1 )  -  10. *X*R(2) / (1 .  +  X*X) 

GO  TO  20 

120  A( 1 )  -  ORD( 1 ) 

B(1)  -  R(1) 

0(1)  -  R(2) 

I- 3 
J  -  K 

II- 2 
12-2 


2k 


125 


130 


140 

150 


DO  130  P  =  12, J 
A(P)  «  R ( I ) 

B(P)  =  R(I+1) 

C(P)  =  R ( 1+2) 

D(P)  =  R(I+3) 

I  -  1+4 
CONTINUE 
J  =  J+1 

A(J)  -  ORD (II) 

11  =  11+1 
B(J)  -  R ( I ) 

IF  (I.EQ.Q)  GO  TO  140 
C(J)  -  R(I+1 ) 

D(J)  =  R(I+2) 

I  =  1+3 
J  -  J+K-1 

12  -  J+2-K 
GO  TO  125 
RETURN 
RETURN  1 
END 


25 


nnnnnfinannTnnnananonnnannnnnnnnonnnoflannnnnnn 


SUBROUTINE  BANDET(M,1DIM,INT,A,N,M1  ,M3,*) 
DIMENSION  M(IDIM,M1) ,INT(IDIM) ,A(1DIM,M3) 


M - AN  NXM1  MATRIX  FOR  STORING  LOWER  TRIAMGUIAR 

MATRIX  OF  LU  DECOMPOSITION  OF  A 

1  NT - AN  NX1  VECTOR  FOR  RECORDING  ROW  INTERCHANGES 

DURING  DECOMPOSITION 

A - AN  NX(M1+M2+1)  MATRIX  WHOSE  COLUMNS  ARE  THE  DIAGONALS 

OF  C,  THE  BAND  MATRIX  BEING  DECOMPOSED 
A(* ,  1 )  -  A(*,M1)  ARE  SUBDIAGONALS  Of  C 
A(*,M1+1)  IS  DIAGONAL  OF  C 

A(*,M1+2)  -  A(* .M14M2+1 )  ARE  SUPERDIAGONALS  OF  C 

N - NUMBER  OF  ROWS  IN  A 

Ml - NUMBER  OF  SUBDIAGONALS  IN  C 

M3 - TOTAL  NUMBER  OF  DIAGONALS  IN  C  ,  I.E.  WIDTH  OF  BAND 

M3  -  Ml  (//  SUBDIAGS)  +  M2  (If  SUPERDIAGS)  +1 

BANDET  AND  BANDSL  ARE  TWO  SUBROUTINES  WHICH  SOLVE  C*X  -  B 
WHEN  C  IS  AN  UNSYMMETRIC  BAND  MATRIX  (  THEY  WILL  WORK  WITH 
SYMMETRIC  BAND  MATRICES  BUT  TAKE  NO  ADVANTAGE  OF  THEIR 
STRUCTURE) . 

C  HAS  Ml  SUBDIAGONALS  AND  M2  SUPERDIAGONALS . 

THE  MATRIX  C  IS  TRANSFORMED  TO  A  BY  MAKING  EACH  DIAGONAL  OF 
C  A  COLUMN  OF  A.  THUS  A  IS  NX(M14M2+1)  WHEN  C  IS  NXN. 

A  TYPICAL  A  IS  PICTURED  BELOW.  C  HERE  IS  4X4,  WITH 

2  SUBDIAGONALS  AND  1  SUPERDIAGONAL 

0  0  C(1 , 1)  C(1 ,2) 

0  C ( 2 . 1 )  C(2,2)  C(2,3) 

C(3,1)  C(3,2)  C  ( 3 , 3)  C(3 ,4) 

C(4,2)  C(4, 3)  C(4 ,4)  0 

THIS  TRANSFORMATION  IS  THE  FOLLOWING,  ASSUMING  THAT  C(I,J) 
IS  A  BAND  ELEMENT  IN  C: 

C(I,J)  — >  A(I,M1+1+(J-I)) 

ALL  OTHER  ELEMENTS  OF  A  ARE  0 

BANDET  FINDS  THE  LU  DECOMPOSITION  OF  A,  STORING 
THE  LOWER  TRIANGULAR  MATRIX  IN  M  AND  NXM1  MATRIX, 

AND  OVERWRITING  THE  UPPER  TRIANGULAR  MATRIX  INTO  A. 

BANDSL  USES  THIS  DECOMPOSITION  TO  SOLVE  A*X  *=  B  WHERE 
THE  RIGHTHAND  SIDE  IS  INPUT  IN  THE  VECTOR  B,  AND  X  IS 
OUTPUT  IN  THE  VECTOR  B. 

THESE  ROUTINES  WERE  TRANSLATED  FROM  THOSE  PRESENTED  BY 
J.  WILKINSON  IN  NUMERISCHE  MATHEMATIK  VOL  9,  P279 
C  TRANSLATOR:  BARBARA  RYDER 
C 
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REAL  M 
25  L-M1 

DO  AO  1-1, Ml 

K2-M 1+2-1 
DO  50  J-K2.M3 
50  «(I,J-L)-A(I,J) 

L-L-1 
K2-M3-L 
DO  60  J-K2.M3 
60  A(I , J)-0. 0 

AO  CONTINUE 
A  5  L-M1 

DO  70  K-1  ,N 
X-A(K, 1 ) 

I-K 

K2-K+1 

IF  (L.LT.N)  L-L+1 

72  IF  (L.LT.K2)  00  TO  81 

79  DO  80  J-K2  ,L 

IF  (ABS(A(J,1))-ABS(X))  80,80,82 

82  X-A(J, 1 ) 

I-J 

80  CONTINUE 

81  INT(K)-I 

IF  (X)  73,75,73 

73  IF  (I-K)  77,78,77 

77  DO  90  J-1.M3 
X«A(K,J) 

A(K, J)-A(I , J) 

90  A(I , J)-X 

78  IF  (L.LT.K2)  GO  TO  70 

83  DO  100  J-K2  ,L 
M(K,J-K)-A(J,1)/A(K,1) 

X-M(K, J-K) 

DO  110  JJ-2.M3 

110  A(J,JJ-1)«A(J,JJ)-A(K,JJ)*X 

100  A(J,M3)-0.0 

70  CONTINUE 
RETURN 

75  RETURN  1 
END 


» 


1 
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SUBROUTINE  BANDSL  (M,IDIM,INT,B,A,N,M1  ,M3) 
DIMENSION  INT(IDIM)  ,A(IDIM,M3)  ,M(IDIM,M1)  ,B(IDIM) 


ALL  PARAMETERS  SAME  AS  IN  BANDET  EXCEPT  FOR: 

B - RIGHTHAND  SIDE  OF  LINEAR  SYSTEM  C*X  -  B 

SOLUTION  IS  RETURNED  IN  B 

REAL  M 
INTEGER  W 
L-M1 

DO  10  K-1,N 
I-INT(K) 

IF  (I-K)  11,12,11 

11  X-B(K) 

B(K)-B(I) 

B(I)-X 

12  K2-K+1 

IF  (L.LT.N)  L-Lfl 

14  IF  (L.LT.K2)  GO  TO  10 

15  DO  20  I-K2.L 
X-M(K.I-K) 

20  B(I)-B(I)-X*B(K) 

10  CONTINUE 
L-1 

DO  30  II-1,M 

I-N+1-II 

X-B(I) 

W-I-1 

IF  (L-1)  32,33,32 

32  DO  40  K»2,L 

40  X*»X-A(I  ,K)*B(K+W) 

33  B(I)“X/A(I, 1 ) 

IF  (L-M3)  31,30,30 
31  L-L+1 

30  CONTINUE 
RETURN 
END 
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3  A  Mehlum*  s  method 


Even  Mehlum  (1969),  in  his  Ph.D.  dissertation,  presented  an  algorithm 
for  computing  an  approximation  to  nonlinear  interpolating  splines  which 
he  has  implemented  in  a  computer  program  called  KlfRGLA  ( from  Norwegian 
"KURveGLAtting"  --  curve  fitting).  KURGLA  is  the  basic  subroutine  in  the 
ship  fairing  program  included  in  the  AUTOKON  software  which  has  been  dev¬ 
eloped  at  the  Central  Institute  for  Industrial  Research  in  Oslo,  Norway. 

The  AUTOKON  computer  software  is  used  for  numerical  control  of  drawing 
machines  and  flame  cutters  in  ship-building.  The  ship  fairing  program 
has  been  in  use  since  1965  and  by  1969  had  been  used  for  fairing  about 
150  hulls. 

Mehlum  starts  with  a  long  derivation  of  the  equation 
2 

|  =  p  sin  (e  -  cp)  +  V  (3*10) 

which  is  yet  another  form  of  Kirchhoff's  kinetic  analogue  (3  A).  He  also 
proves  the  following: 

Theorem  (3-3):  Between  neighboring  knots  of  a  nonlinear  spline,  there  is 
always  a  direction  along  which  the  curvature  varies  linearly. 

Proof:  This  follows  from  the  fact  that  curvature  is  proportional 
to  the  bending  moment  which  can  be  expressed  as  a  linear  function 
of  x  and  y  .  (See  Equation  (2.6)). 

Mehlum  makes  the  assumption  that  m  =  0  ,  so  that  the  direction 
along  which  the  curvature  varies  linearly  coincides  with  the  x-axis.  The 
curvature  is  then  represented  as  a  piecewise-linear  continuous  function 
which  changes  slope  only  at  the  knots.  Mehlum  further  simplifies  the 
problem  by  representing  the  linear  curvature  between  each  knot  by  a  series 
of  steps  giving  a  piecewise  constant  function  (noncontinuous,  of  course) 


lor  the  curvature.  In  other  words,  he  approximates  ',h<*  solution  to  a 
somewhat  different  problem  tnan  ( *’ .  10)  by  a  series  ol'  arcs  oi  circles 
whi  cm  iorm  a  continuous  1  tine  Lion  with  a  continuous  first  derivative,  but 
with  discontinuous  second  derivatives.  Combining  this  representation 
with  ('’.10)  and  integrating  yields  two  constants  of  integration  to  be 
determined  within  each  panel  lor  a  given  set  of  curvatures  specified  at 
each  knot.  His  algorithm  thus  breaks  into  an  inner  and  an  outer  iteration. 
The  outer  iteration  first  guesses  a  set  of  knot  curvatures  and  then  tries 
to  improve  them  based  upon  how  close  the  resulting  curves  fit  the  data 
and  specified  ena  slopes.  (He  is  not  trying  to  produce  natural 
splines.)  The  inner  iteration  then  uses  the  current  values  for  the  knot 
curvatures  and  determines  the  constants  of  integration  for  each  panel 
through  some  simple  recurrence  formulas  so  that  the  function  interpolates 
the  end  knots.  The  output  of  Mehlum’s  subroutine  is  in  the  form  of  center 
coordinates  and  radii  for  the  circles  thus  obtained. 

Mehlum  claims  that  many  numerical  control  machines  can  use  output  of 
this  form  directly. 
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4  .  The  Discrete  Natural  Cubic  Interpolating  Spline 

The  data  points  =  Y(X^)  ,  i  =  1  ,  ...  ,  n  (vhere 

X^  <  ...  <  Xn)  ,  may  be  interpolated  as  follows:  Let 


H.  =  X,  -  X.  ,  i  =  1,  ...  ,  n-1  . 
i  i+1  i 


Assume  that  each  ol  the  is  an  integral  multiple  of  the  mesh  size  parn. 

meter  h  .  The  interval  [X^  -  h  ,  X^  +  h]  is  divided  into  a  partition 

x0  -  *1  >  •••  >  xm  ’  xm+l  >  where 


xi  =  X^  +  (i-l)h  ,  i  =  0  ,  ...  ,  m  +  1  ,  and 

m  =  (xn  -  X1)/h  +  1  . 

Denote  by  X  the  set  ol  abscissa  data  X^  ,  i  =  1  ,  ...  ,  n  ,  and 
let 

\  =  0  , 
i-1 

Ln-  =  H.  ,  i  =  2  ,  ...  ,  n  . 

0=1 

Tin  discrete  natural  cubic  spline  interpolating  the  points  (X^  ,  Y^)  , 
i  =  1  ,  ...  ,  n  i r  defined  to  be  tne  set  of  points  (x^  ,  y^)  , 
i  =  0  ,  ...  ,  m  +  1  ,  which  minimize  the  value  of  S  where 

>2 

h 
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subject  to  the  constraints  of  interpolation, 

yt  =  Yj  ,  if  i  =  Lj/h  +  1  . 

Let  f>  denote  the  central  difference  operator.  The  raesh  ordinates 
y^  ,  i  -  2  ,  ...  ,  m-1  ,  are  computed  by  solving  the  linear  system  of 
equations 

ft^y^  =  0  if  x^X  and  0  <  i  <  m  ,  (4.1) 

and 

.2  ,2  ~ 

6  yi  -  6  ym  •  0  , 

which  arise  by  setting  ds/dy^  =  0  ,  i=0,2, . . . ,m+l  (excluding  every 

i  =  Lj/h  ,  for  j=l,...,n),  subject  to  the  conditions 

y^^  =  Yj  if  i  =  Lj/h  ,  i  =  1, . . . ,m  . 

The  values  of  y~  and  y  , ,  are  not  explicitly  computed,  but  are 
u  m+i 

introduced  into  the  formulation  to  accommodate  the  conditions  on  the 
second  differences  f  t  the  end  nodes  and  the  first  and  last  of  equations 
(4.1).  Notice  the  similarity  between  this  formulation  and  that  of  the 
continuous  natural  cubic  spline. 

A  useful  representation  of  the  above  linear  system  is  given  by 

Ay  =  b  , 

where  y  =  (y,  ,  ...  ,  y  )  ,  the  matrix  A  has  the  structure 
~  x  in 
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10  0 
5-4  1 

1-4  6-1+  l 
l-i+  6-4  1 


1-4  6-4  1 

0  0  10  0 
1-4  6-4  1 


1-4  6-4  1 

0  0  10  0 
1-4  6-4  1 

•  •  •  •  • 

•  •  •  •  • 

1-4  6-4  1 
1-4  5-2 

0  0  1 


and  b  =  (Yx  ,  0  ,  ,0,^,0,  ...  ,  0  ,  y  .  ...  0  ,  Yn)T  .  The 

row. j  of  A  with  onr  a  on  the  main  diagonal  cori\  spond  to  thr  nonzero 
dementi:  in  the  vector  b  .  The  second  row  of  A  combines  the  two 


1 

j 

A 

i 


i 

S 


i 


j 

i 

i 


*  ouation, 


The  (m-l)th  row  of  A  is  formed  similarly. 


The  following  analysis  provides  a  way  of  solving  Ay  =  b  by  solving 
a  symmetric  band  system  of  smaller  order  than  A  .  A  bound  on  the  con¬ 
dition  number  of  this  symmetric  matrix  is  determined  and  found  to  be 
independent  of  the  order  of  the  matrix  (total  number  of  mesh  points  and 
knots)  but  highly  dependent  upon  the  number  of  mesh  points  between  knots. 

A  symmetric  matrix  closely  related  to  the  matrix  A  can  be  deter¬ 
mined  as  follows.  Let  the  matrix  B  have  the  same  structure  as  A 
except  that  the  columns  of  B  having  ones  on  the  main  diagonal  have 
all  of  the  other  components  set  to  zero.  Thus,  the  vector  y  can  be 
computed  by  solving  the  system 

By  =  b ' 

where  the  vector  b  '  is  an  appropriate  modification  of  the  vector  b  . 

Now  since  B  is  a  symmetric  matrix,  its  condition  number  is  the  maximum 
ratio  of  the  magnitudes  of  its  eigenvalues.  By  expanding  the  determinant 
of  B  -  XI  along  any  row  or  column  of  B  having  a  one  on  the  main 
diagonal,  it  is  obvious  that  X  =  1  is  a  multiple  eigenvalue  of  B  . 

These  eigenvalues  may  be  eliminated  from  B  by  dropping  the  rows  and 
columns  with  1  on  the  main  diagonal,  leaving  the  matrix 


where  each  eigenvalue  of  B*  is  an  eigenvalue  of  B  ,  In  practice,  the 
unknown  components  of  the  vector  y  are  computed  by  solving  a  linear 
system  having  B*  as  the  coefficient  matrix.  Now 

B*  =  C  +  D  ,  where 
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and  each  is  a  square  matrix  of  order  JL/h-1  of  the  form 


The  matrix  D  has  eigenvalues  0  and  2  .  If  the  eigenvalues 


X.(B*)  and  \i(C)  ,  i  =  1 


?  •  •  •  5 


m  -  n 


are  arranged  in  a  non- 


increasing  order,  then  by  a  corollary  ol‘  the  Courant -Fisher  theorem 
(see  Wilkinson  (1965 ),  p.  101), 

X  (B*)  >  X^C)  ,  i  =  1  ,  ,  m  -  n  . 


The  eigenvalues  of  C  are  simply  those  of  the  Ci  ,  i  =  1  ,  ...  ,  m  . 

2 

Now  the  Ci  can  be  written  as  ,  where 


The  characteristic  polynomials  of  the  upper-left  hand  minors  of 
satisfy  the  same  three-term  recurrence  relation  as  the  Chebyshev  poly¬ 
nomials,  namely 


T.^X)  =  0  , 

T0(X)  =  1  , 

Tk+i(X)  =  (-2-X)Tk(X)  -  T^X)  ,  k  =  0,1,...  . 
It  follows  that 


where 


Xj(D.)  =  -2 ( 1-cos  ^  )  ,  j  =  1 


v  =  H^/h-1  .  Thus, 

min  Xj(C.)  =  4(1  -  cos  ^  f  , 

J 


v  , 


and  therefore, 
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min  |X.(B)|  >  min  4(1  -  cos  —  ) 
i  j  i  Hi 


By  Gerschgorin's  theorem, 


max  I  X  (B) |  <  l6  . 


This  proves 


Theorem  4.1: 


cond  (B)  < 


min  (1  -  cos  ~  )c 
i  Hi 


_ *  ..  <  *«,  i 

min  (  jj-  )  i  \ 

i  i 


Thus,  for  many  problems,  the  linear  system  resulting  from  the  finite 
difference  formulation  is  reasonably  well-conditioned.  More  importantly, 
the  bound  on  the  condition  number  of  B  is  independent  of  n  ,  the 
number  of  data  points. 

The  purpose  of  the  following  discussion  is  to  analyze  the  difference 
between  y  and  the  restriction  of  the  continuous  natural  cubic  spline 
function  interpolating  the  points  (Xi  ,  Y^,  i  =  1  ,  ,  n  ,  to 

the  abscissae  xi  ,  i  =  ],...,  m  .  We  call  this  difference  the  dis¬ 
cretization  error  of  y  and  will  now  bound  a  certain  norm  of  it. 

Consider  the  neighborhood  of  : 


.  .  -2  -1 


1  2 


w 


For  convenience,  the  mesh  nodes  are  now  numbered  ...  ,  -2  ,  -1  ,  0  ,  1  , 

2  ,  ,  such  that  Xq  =  . 

k 

Since  6  y^  =  0  for  i  =  . . .  ,  -2  ,  -1  ,  it  follows  that  the 
points  ,  for  i  =  ...  ,  -2  ,  -1  ,  0  ,  1  lie  on  a  cubic  polynomial. 

Similarly,  the  points  y^  ,  for  i  =  -1,0,  1,2,  ...  ,  lie  on  a 

(generally  different)  cubic  polynomial.  These  two  cubics  coincide  at 
the  points  -1,0,  and  1  .  Let 

^i  =  ai  +  ^i*  +  ci^  +  d^  ,  where  ?  =  x  -  X^  , 

denote  such  a  cubic  for  x«  [X.^  -  h  ,  X^+1  +  h]  .  The  coefficients  a.^  , 

bi  ,  ci  and  d.^  ,  i  =  l,...  ,  n  -  1  ,  can  actually  be  computed  by 

observing  the  following:  ty."  is  linear.  Let  *^(x)  =  ty^(x)  if 
xe[X^,X^+^]  ,  i  —  1  ,  ...  ,n  —  1.  Thus, 

V(x)  *  V  *  *  H71"  (Ii+i  - 

where  fL  was  previously  defined  as  Hi  =  Xi+1  -  Xi  ,  i=l  ,  . . .  ,  n-1  . 
Note  that  (x)  is  continuous  at  each  knot  since 


for  each  cubic  ^  and  in  the  interval  [X  -  h  ,  X^  +  h] 

Integrating  (4.2)  gives  2 

V'(x)  =  \  *  y.'tx  -  h)  ♦  (y"+1  -  y”)  ■  '*  ‘  V 


2H. 


where  +Y^  is  a  constant  of  integration,  and 


?4(x)  =  Y±  +  +Yt'(x  -  X±)  +  Y" 


(x  -  X,)‘ 


y  r 


(x  -  xi)3  • 


(4.3) 


(4.4) 
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Here  +Y^  denotes  ^(Xi+)  .  Since  ty'  is  not  necessarily  continuous 
at  the  knots,  will  be  used  to  denote  *^/(Xi~)  .  Evaluating  (4.4) 


st  Xi+1  gives 


V  Yi+1  -  Y1 
i  "  H, 


i  3 


(4.5) 


Replacing  i  by  i  -  1  in  (4.3)  and  evaluating  it  at  Xi  gives 


=  Vi.!  +  d'i.l  +Y'i> 


(■>.6) 


Now,  for  y  , 


fiy0  y,  -  y  ,  2 

-ST-  =  -5T"  =  bi  +  dih 


-  V  +  dih2  , 


and  for  y  , 


k  [ai-i +  bi-i  <Hi.i +  h> +  °i-x  (Hi-i +  h>2  +  ai-i  (Hi.i +  h>' 

ai-l  -  bi-l  <Hi-l  -  b>  -  ci<Hi-l  -  b>2  -  V*l-1  -  b)3J 


=  b.  .  +  2c.  H.  .  +  3d.  _  H?  .  +  d.  ,  hc 
l-l  i-l  l-l  l-l  l-l  l-l 


=  V  +  di-ih  • 
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* 


» 


Thus, 

"V  =  +V  +  h2(di  “  d^l)  •  (4-7) 

Combining  (4.7)  with  (4.5)  and  (4.6)  gives 


^'i-l  Hi-l  +  2Y"i  (Hi-l  +  Hi)  +  Y"i+1  Hi  +  6h2(di  "  di-l)  (4‘8) 
=  6(A^+^  “  ^i^  >  i=2,...,n-l, 

where 

-  (y±  -  Y^_^)/H^_^  »  i  =  1  ,  ...  ,  n  . 

2 

The  conditions  6  y.  =0  at  the  ends  require  that  Y"^  =  Y"n  =  0  . 
Equations  (4.4)  and  (4.5)  give 


ai 


Yi  * 


b.  = 


Vl  -  Yi 

Hi 


c,  =  rj2  ,  and 


d.  = 


i+1 


-  Y". 


6h. 


y  1  —  1  J  •••  ,  n  -  1  . 


(4.9) 

(4.10) 

(4.11) 

(4.12) 


Substituting  (4.12)  into  (4.8)  gives  a  linear  system  of  equations  in  the 
unltnowns  Y"  ±  ,  i  =  2  ,  . . .  ,  n  -  1  ,  the  solution  of  which  can  be  used 
to  compute  the  coefficients  ,  b^  ,  c^  ,  and  d^  ,  i  =  1  ,  ...  >  n  -  1  > 
using  (4.9)  -  (4.12)  . 


» 

4l 


t 


t 


A  practical  method  for  computing  the  coefficients  of  the  continuous 
cubic  spline  interpolating  {X^  ,  Y^)  ,  i  =  0  ,  . . .  ,  n  ,  can  be  formulated 
by  following  the  above  steps  and  replacing  (4.6)  with  the  continuity 
conditions 

-Y'  -  \  ,  i  =  1  ,  ,  n  -  1  . 

This  leads  to  a  system  of  equations  for  Y*  which  is  identical  to  (4.8) 

2 

without  the  6h  (d^  -  d^)  te:m.  Thus  we  have  the  same  system  of  equations 

n  2 

for  Y^  in  the  discrete  case  as  in  the  continuous  case,  except  for  an  0(h  ) 
perturbation  of  the  coefficient  matrix. 

Let  s"  denote  the  vector  of  second  derivatives  (Yp  for  the 
continuous  case  and  let  y"  denote  the  vector  of  second  derivatives  for 

Also,  let 


the  discrete  case. 


" 

at  =  y  -  s 


2(H1+H2)  H2 

h2  ^ 

Hj  2(H5+H4) 


H 


n-2 


H  0 
n-2 

2(H  _+H  J 
n-2  n-1' 


and 
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Let  ^  =  6Ui+1  ~  A^)  .  Then  from  (4.8), 
(A  +  6A)  (s/7  +  a)  =  8  . 


a  =  [(A  +  fiA)"1  -  A-1]  0  . 

Setting  B  =  A  +  6A  in  the  identity 
B'1  -  A-1  =  A_1(A  -  B)  B"1  , 
and  substituting  gives 

a  =  -A_16A(A  +  5A)"1  B  =  -A_1ftA(s'+  a)  . 
Taking  spectral  norms  gives 
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11*11  <  IIA’II  ||6A|f  ||s'  +  *| 


II  all  ,  %  m\ 

■  <  cond(A)  - 

II  y"ll  "  II  a  || 

Now,  by  Gershgorin's  theorem, 


Thus, 


min  (H^  +  Hi+1)  <  X(A)  <  3  max  (H^  +  Hi+1) 
i  ”  i 

X(fiA)  <  max  2h2  +  rp —  J  <  max  . 

”  i  \  1  i+l  /  ”  i  l 

cond(A)  <  3  max  (H.  +  Hi+1)/min  (H.  +  H.+1)  ,  and 
i  i 


||y//j|  ~  (min  IL)  min(HL 

Since  the  uniform  norm  of  a  vector  is  at  most  equal  to  the  spectral  norm, 


4h2||y"|| 

-  H  .  min(H?  +  H..J 
0  mm  .  '  l  i+l' 

l 


,  J=l,.,,  ,n->l. 


The  coefficients  a.  .  b.  ,  c.  ,  and  d.  are  linear  functions  of  the 

,  hence  the  difference  between  the  continuous  and  discrete  splines 
2 

is  0(h  )  .  More  precisely,  denoting  the  continuous  spline  by  s(x)  , 
equations  (4.9)  -  (4.12)  yield 


Theorem  4.2: 


,  <wVuh2 

max  |s(x)  -7W|  <T  BTT) 

x  ^min  .  '  l  i+l' 

i 


,  x6  [XQ,Xn]  . 
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Theorem  4 .2  tells  us  that  the  discrete  natural  cubic  interpolating 

spline  coincides  with  the  continuous  natural  cubic  spline  in  the  limit 

as  the  mesh  size  h  is  decreased  to  zero.  As  the  mesh  size  is  decreased, 

the  error  in  approximating  the  continuous  spline  decreases  at  least  as 

fast  as  the  square  of  the  mesh  size.  However,  Theorem  4.1  tells  us  that 

the  condition  number  of  the  coefficient  matrix  of  the  linear  system  may 

4 

increase  as  fast  as  h  .  On  a  given  machine,  it  should  thus  be  possible 
to  calculate  a  reasonable  discrete  spline  approximation;  however,  problems 
with  fine  meshs  could  be  difficult,  if  not  impossible,  to  solve. 

The  work  required  to  compute  discrete  cubic  splines  is  consider¬ 
ably  more  than  that  required  to  compute  continuous  cubic  splines.  The 
primary  reason  for  studying  them  here  is  to  gain  insight  into  the  algo¬ 
rithm  for  computing  discrete  approximations  to  nonlinear  splines,  dis¬ 
cussed  in  the  following  section.  In  view  of  the  difficulty  in  analyzing 
linear  discrete  splines,  it  comes  as  no  surprise  that  the  nonlinear 
case  has  not  been  analyzed.  The  proofs  of  this  section  rely  upon  special 
properties  of  both  the  discrete  and  continuous  cubic  splines.  Hence, 
few,  if  any,  of  the  arguments  would  carry  over  to  the  nonlinear  case. 

Discrete  cubic  splines  and  generalized  discrete  splines  have  been 
characterized  by  Mangasarian  and  Schumaker  (l$Tfl).  They  prove  existence 
and  uniqueness.  Some  weak  convergence  results  for  discrete  splines  have 
been  published  by  Daniel  (lS^l).  Theorem  4.2  appears  to  be  the  only 
rate-of-convergence  result  for  discrete  splines. 
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In  this  section,  an  iterative  method  is  presented  for  computing 


a  discrete  approximation  to  the  nonlinear  spline  interpolating  the  data 
points  (Xi,Yi)  ,  i=l,...,n  . 

As  in  the  previous  section,  we  assume  that  <  X2  <  . . .  <  Xq  , 
and  the  interval  [X^  -  h  ,  Xn  +  h]  is  partitioned  into  a  uniform  mesh 


x  ,,  ,  where 
m+1  ’ 


xi  =  X1  +  (i-l)h  ,  i=0,l, . . . ,m+l  , 

and 


m  =  (xn  -  Xx)/h  +  1  , 


for  some  mesh  size  parameter  h  .  For  simplicity,  assume  that  the 
X±  (i=l, . . .  ,n)  are  equally  spaced  with  k  mesh  intervals  between  con¬ 
secutive  data  points,  i.e., 

k  =  (X.+1  -  Xi)/h  ,  i=l, . . . ,n-l  . 

We  compute  a  discrete  approximation  y  =  [y^,...,y  ]  *'unc" 

tion  y  which  minimizes  the  functional 


The  functional  (5.1)  is  approximated  by  substituting  a  summation  for 
the  integral, 

h"262y.  «  y//(xi)  , 

and 

(yi+i  -  «  y/(xi)  . 
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This  gives  the  function 


,  (y) .  f1  (y^i  -  gyi : 

{ri  11  *  (yi.i  '  yi-i )2/l,h2)5 


to  be  minimized  subject  to  the  interpolation  constraints 


(5.2) 


yi  =  Yj  ’  if  Xi  Xj  ’  ‘ 

Again,  the  fictitious  ordinates  yQ  and  ym+1  are  introduced  to  accom 
modate  the  first  and  last  terms  of  (5.2). 

A  necessary  condition  for  E^(y)  to  achieve  a  minimum  is 
§y~  =  ^  *  for  a^-  i  /  1  (mod  k)  . 

This  yields  the  nonlinear  system  of  equations 
t>\  =  0  , 

- 0  - 

'i-i(yi '  2yi-i  *  yi-2>  -  e?i(yi+i  -  2yt  +  (5.3) 

+  !M(yit2  -  2yi+l  +  yi> 

•Xl-l(yi  -  yi-2>  *  Wyi*2  -  yi>  -  0  . 

where 

<i  =  f1  +  <y1+1  -  , 

h  ’  (yi+i  -  2yi  +  yi-i>2[i  *  <yi+1  -  yt ./Ah8]'772  , 

(for  i  =  2, . . .  ,m-l  ,  except  i  =  1  (mod  k))  . 

The  System  (5*3)  can  be  solved  in  many  ways.  A  simple,  though 
probably  not  the  fastest  way  is  a  Picard-type  iteration  described  as 


kl 


This  algorithm  is  implemented  in  the  Fortran  subroutine  SPLINE 
found  at  the  end  of  this  section. 

A  Cholesky  method  is  used  for  the  solution  of  (5.4).  This  choice 
is  based  on  computational  experience.  The  coefficient  matrix  of  (5.4)  is 
usually  positive  definite.  Cases  in  which  the  coefficient  matrix  is  in¬ 
definite  appear  to  correspond  to  those  in  which  either  no  single-valued 
solution  exists  or  the  outer  iteration  becomes  unstable  due  to  large 
slopes  in  the  solution. 
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Various  experiments  have  been  curried  out  using  both  SPLINE  and 
MEC  given  in  Section  3.3.  The  first  of  these  used  the  Woodford  data 
given  in  Section  3.3.  For  this  case,  n  •=  7  •  Figure  5-1  shows  times 
required  to  compute  splines  of  similar  accuracy  on  an  IBM  3^0/ 67  for 
various  values  of  k  .  These  runs  used  EPS  values  of  10~^  .  Figure 
5.1  indicates  that  both  SPLINE  and  MEC  require  time  proportional  to  k  , 
though  spline  appears  to  be  about  an  order  of  magnitude  faster  for  this 
problem.  For  values  of  k  greater  than  40,  MEC  failed  -  though  this 
could  be  partly  due  to  the  fact  that  MEC  is  programmed  in  single  precision. 
SPLINE  worked  for  values  of  k  as  large  as  140.  Values  of  (5.2)  are 
given  in  Table  5.1  for  solutions  given  by  MEC,  SPLINE  and  the  natural 
continuous  cubic  spline,  for  Woodford's  data  and  various  values  of  k  . 

The  values  of  Eh(y)  for  the  discrete  cubic  spline  agreed  with  those 
for  the  continuous  cubic  spline  to  three  significant  digits.  The  last 
column  of  Tab’e  5-1  indicates  that  the  results  of  SPLINE  and  MEC  con¬ 
verge  uniformly  to  one  another  as  k  increases.  The  nonlinear  spline 
interpolating  Woodford's  data  is  shown  in  Figure  5.2  as  plotted  by  a 
Versatec  electrostatic  plotter. 

An  experiment  to  determine  the  dependency  of  computation  time 
upon  n  ,  the  number  of  data  points,  was  constructed  as  follows.  The 
number  of  mesh  intervals  between  adjacent  data  points  was  kept  constant 
at  k  =  10  .  The  mesh  size  was  h  =  .1  giving  -  x^  =  kh  =  1  , 

i  =  l,...,n-l  .  The  ordinates  were  chosen  as 
\\  =  (i  mod  2)/5  ,  i  =  l,...,n  , 

to  insure  a  solution  having  relatively  small  slopes  and  curvatures. 

The  parameters  EPS  were  set  to  10~^  .  Measured  CPU  times  are  shown 
in  Figure  5.3  for  both  SPLINE  and  MEC.  The  two  values  of  n  shown  for 


49 


MEC  in  Figure  5 .  j  are  n  =  5  and  n  =  10  .  From  these  two  times,  MEC 
does  not  appear  to  r  quire  time  proportional  to  n  .  For  values  of  n 
greater  than  10  ,  MEC  failed  to  converge.  SPLINE,  on  the  other  hand, 
required  relatively  small  amounts  of  time  which  are  proportional  to  n  . 
SPLINE  performed  well  for  values  of  n  as  large  as  100.  The  only  res¬ 
trictions  on  n  for  SPLINE  appear  to  be  available  core  storage  and 
computation  time. 
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10 

2.55 

2.52 

2  .69 

.020 

^0 

2.55 

2.55 

2.69 

.011 

50 

2.55 

2.55 

2.70 

.0070 

to 

2.55 

2.55 

2.70 

1 _ 

.0051 

Table  5.1.  Values  of  E^(y)  ^or  Woodford's  data  for 
MEC,  SPLINE,  and  the  cubic  spline,  and  a  uniform-norm  comparison 

of  MEC  with  SPLINE. 


Figure  5.2.  Nonlinear  spline  interpolating  Woodford's  data 
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SUBROUTINE  SPLINE  (N,  ORD,  K,  H,  Y,  EPS,  *) 

INTEGER  N,  K 

DOUBLE  PRECISION  ORD(N),  H,  Y(l),  EPS 
C 

C  THIS  SUBROUTINE  COMPUTES  THE  NONLINEAR  INTERPOLATING  SPLINE  WHICH 
C  PASSES  THROUGH  THE  N  DATA  POINTS  WITH  ORDINATES  ORD  (I), 

C  (1-1,...  ,N),  AND  EQUIDISTANT  ABSCISSAE  X(I),  (1-1 . N) ,  SATISFYING 

C  X(1).LT.X(2).LT . LT.X(N)  ,  AND  X(I+1)-X(I)  -  K*H,  (1-1 , . . .  ,N-1 )  . 

C  THE  DISCRETE  APPROXIMATION  Y(I),  (1-1 , . . . , (N-1 )*K+1)  IS  RETURNED 

C  WITH  ORDINATES  CORRESPONDING  TO  THE  MESH  ABSCISSA  X(1)+(I-1)*H,  WHERE 
C  H  IS  THE  SPECIFIED  MESH  SIZE  PARAMETER.  K  IS  THE  NUMBER  OF  MESH 
C  INTERVALS  BETWEEN  SUCCESSIVE  DATA  POINTS.  THE  INPUT  PARAMETER  EPS 
C  IS  USED  TO  DETERMINE  WHEN  TO  STOP  THE  ITERATIONS.  ROUGHLY,  IF  P 
C  SIGNIFICANT  DIGITS  ARE  DESIRED  IN  A  PROBLEM  WHOSE  SOLUTION  IS  0(1), 
C  THEN  EPS  SHOULD  BE  ABOUT  1.D-P. 

C  THE  ERROR  RETURN  IS  TAKEN  IN  THE  EVENT  THAT  A  RESULTING  LINEAR 

C  SYSTEM  IS  TOO  ILL-CONDITIONED  TO  SOLVE. 

C  K  MUST  BE  AT  LEAST  5. 

C 

C 

DOUBLE  PRECISION  KZ(IOOO),  C(1000,3),  B(1000),  NORM,  H2,  DIF 
DOUBLE  PRECISION  F8H2,  LZ(IOOO),  Y1 
DOUBLE  PRECISION  DABS 
INTEGER  M,  Ml,  MM,  K2,  N1,  I,  J,  L,  J2 
C 

H2  -  H+H 

F8H2  -  5./(8.*H**2) 

M  -  K* (N-1 )+1 
Ml  -  M-1 
K2  -  K-2 
N1  -  N-1 


INITIALIZE  Y 

DO  5  1-1  ,M 
Y(I)  -  0. 

5  CONTINUE 

SET  DATA  ORDINATES  IN  Y 
J  -  1 

DO  10  1-1 ,M,K 
Y(I)  -  ORD(J) 

J  »  J+1 
10  CONTINUE 

INITIALIZE  KZ  TO  ONES  FOR  THE  FIRST  ITERATION 


nnn  o  o  o  o  ooo  ooo 


DO  20  1-2, Ml 
KZ(I)  -  1. 

LZ(I)  -  0. 

20  CONTINUE 

KZ ( 1 ) , KZ (M) , LZ ( 1 ) ,  AND  LZ(M)  ARE  SET  TO  0.  TO  SIMPLIFY  BOOKKEEPING 

KZ(1)  -  0. 

KZ(M)  -  0. 

LZ(1)  -  0. 

LZ (M)  -  0. 

GO  DO  FIRST  ITERATION 
GO  TO  50 

BEGINNING  OF  MAIN  LOOP, 

RECOMPUTE  NEW  KZ  AND  LZ  VECTORS 

30  DO  40  1-2, Ml 

Y1  -  1 .  +  ( (Y (1+1 )-Y (1-1 ) ) /H2)**2 
KZ(I)  -  Y1**(-2.5) 

LZ(I)  -  F8H2*(Y(I+1)-2.*Y(I)+Y(I-1))**2*KZ(I)/Y1 
40  CONTINUE 

SET  UP  COEFFICIENT  MATRIX  AND  RIGHT-HAND  SIDE 

50  L  -  0 

DO  80  1-1 ,N1 
L  -  L  +  1 
MM  -  (I- 1 ) *K 
C(L, 1 )  -  0. 

C(L,2)  -  KZ (MM+1 )  +  LZ(MMfl) 

C(L,  3)  -  KZ  (MM+1 )  +  4.*KZ(MMf2)  +  KZ(MM+3)  -  LZ(MM+1)-LZ(MM+3) 
B(L)  -  2.*(KZ(MM+2)+KZ(MM+1))*Y(MMf1) 

L  -  L+1 
C (L, 1 )  -  0. 

C (L, 2)  -  -2.*(KZ(MM+2)  +  KZ(MM+3)) 

C(L,3)  -  KZ (MM+2)  +  4.*KZ(MMf3)  +  KZ(MMf4)  -  LZ (MMf 2 ) -LZ (MM+4 ) 
B(L)  -  -(KZ (MMf2)+LZ(MM+2)  )*Y (MMf  1) 

L  -  L+1 

IF  (K2.LT.4)  GO  TO  70 
DO  60  J  -  4,K2 

MM  -  (1-1 )*K  +  J 
C (L, 1 )  -  KZ (MM- 1 )  +  LZ(MM-I) 

C(L, 2)  -  -2.*(KZ(MM-1)  +  KZ(MM) ) 

C(L,3)  -  KZ  (MM- 1 )  +  4 .  *KZ  (MM)  +  KZ(MM+1)  -  LZ(MM-1  )-LZ(MMf  1 ) 
B(L)  -  0. 

L  -  L+1 

60  CONTINUE 
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70 


MM  -  l*K 

C  (L,  1 )  -  1C’ (MM- 2)  +  LZ(MM-2) 

C(L,2)  -  -2,*(KZ(MM-1)  +  K2'(MM-2)) 

C(L,3)  -  KZ (MM- 2)  +  4.*KZ(MM-1)  +  KZ(MM)  -  LZ(MM-2)-LZ(MM) 
B(L)  -  -  (KZ  (MM)+LZ  (MM) )  *  Y(MMfl) 

L  -  L+1 

C(L,1)  -  KZ(MM-I)  +  LZ(MM-I) 

C(L, 2)  -  -2.*(KZ(MM-1)  +  KZ (MM) ) 

C(L,3)  -  KZ  (MM- 1 )  +  4  .*KZ  (MM)  +  KZ(MM+1)  -  LZ(MM-1)-LZ(MM+1) 
B(L)  -  2 .  *  (KZ  (MM)+KZ  (MM+1 )  )*Y  (MMf  1 ) 

80  CONTINUE 

SOLVE  5-BAND  LINEAR  SYSTEM 

CALL  CCOMP(L,  3,  1000,  C,  C,  &100) 

CALL  COLVE (L,  3,  1000,  C,  B,  B) 

UPDATE  SOLI. "ION  VECTOR  Y  AND  CHECK  FOR  CONVERGENCE 

L  -  1 
NORM  -  0. 

DO  90  1-1, N1 
DO  90  J-2.K 

J2  -  (1-1 )*K+J 

DIF  -  DABS(Y(J2)  -  B(L) ) 

Y(J2)  -  B(L) 

IF  (DIF. GT. NORM)  NORM  -  DIF 
L  -  L+1 
90  CONTINUE 

IF  (NORM.GT.EPS)  GO  TO  30 
RETURN 
100  RETURN  1 
END 
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SUBROUTINE  CCOMP(N,  M,  IDIM,  C,  L,  *) 

INTEGER  N,  M,  IDIM 

DOUBLE  PRECISION  C(IDIM,  M) ,  L(IDIM,  M) 

C 

C  THIS  SUBROUTINE  FINDS  THE  CHOLESKY  DECOMPOSITION  OF  THE  (2M-1)-BAND 
C  MATRIX  A  OF  ORDER  N  WHERE  THE  MAIN  AND  M-1  SUB  DIAGONALS  OF  A 
C  ARE  STORED  IN  C  AS  SHOWN  (M-3  IS  USED  FOR  THIS  EXAMPLE) 

C 

C  X  X  A(1 , 1 ) 

C  X  A(2, 1)  A(2,2) 

C  A(3, 1)  A(3,2)  A(3, 3) 

C  A(4,2)  A(4 ,3)  A(4,4) 

C  •  «  i 

C  •  .  i 

C 

c 

C  THE  LOWER  TRIANGULAR  PART  OF  THE  CHOLESKY  DECOMPOSITION  A  «  LU, 

C  WHERE  U  IS  L  TRANSPOSE,  IS  RETURNED  IN  L.  THE  RETURN  1  IS  TAKEN 
C  IN  THE  EVENT  THAT  THE  MATRIX  A  IS  FOUND  TO  BE  SINGULAR,  OR 
C  NON  POSITIVE  DEFINITE.  C  AND  L  MAY  BE  THE  SAME  ARRAYS, 

C  IF  DESIRED. 

C  THIS  SUBROUTINE  IS  A  PARTIAL  TRANSLATION  OF  THE  ALGOL  60 

C  PROCEDURE  CHOBANDDET  BY  R.S.  MARTIN  AND  J.H.  WILKINSON,  FOUND 
C  IN  THE  HANDBOOK  FOR  AUTOMATIC  COMPUTATION,  VOLUME  II  (LINEAR 
C  ALGEBRA),  1971,  J.H.  WILKINSON  AND  C.  REINSCH  (EDITORS). 

C 

C 

INTEGER  I,  P,  R,  S,  Q 
DOUBLE  PRECISION  Y 
DOUBLE  PRECISION  DSQRT 
C 

DO  50  1=1, N 

P  =  MAX0(1,M-I+1) 

R  =  I-M+P 
DO  40  J=P,M 
S  =  J-1 
Q  =  M-J+P 
Y  *  C(I,J) 

IF  (P.GT.S)  GO  TO  20 
DO  10  K=P, S 

Y  =  Y  -  L(I,K)*L(R,Q) 

Q  -  0+1 

10  CONTINUE 

20  IF  (J.NE.M)  GO  TO  30 

IF  (Y.LE.O.)  RETURN  1 
1.(1 ,  J)  =  1. /DSQRT  (Y) 

GO  TO  40 

30  L(I , J)  =  Y*L(R,M) 

R  -  R+1 
40  CONTINUE 
50  CONTINUE 
RETURN 
END 
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SUBROUTINE  COLVE  (N,  M,  IDIM,  L,  B,  X) 

INTEGER  N,  M,  IDIM 

DOUBLE  PRECISION  L(IDIM,  M) ,  B(1),  X(1) 

C 

C  THIS  SUBROUTINE  SOLVES  THE  (2M-1)-BAND  LINEAR  SYSTEM  OF 
C  EQUATIONS 
C 

C  AX  -  B 

C 

C  USING  THE  LOWER  PART  OF  THE  CHOLESKY  FACTORIZATION  L  OF  A, 

C  WHICH  HAS  BEEN  COMPUTED  BY  CCOMP. 

C  B  AND  X  MAY  BE  THE  SAME  ARRAYS,  IF  DESIRED. 

C  THIS  SUBROUTINE  IS  A  PARTIAL  TRANSLATION  OF  THE  ALGOL  60 

C  PROCEDURE  CHOBANDSOL  BY  R.S.  MARTIN  AND  J.H.  WILKINSON,  FOUND 
C  IN  THE  HANDBOOK  FOR  AUTOMATIC  COMPUTATION,  VOLUME  II  (LINEAR 
C  ALGEBRA),  1971,  J.H.  WILKINSON  AND  C.  REINSCH  (EDITORS). 

C 

C 

INTEGER  S,  J,  P,  Q,  NP1 ,  K1,  II 
DOUBLE  PRECISION  Y 

SOLUTION  OF  LY  -  B 

DO  20  1-1, N 

P  -  MIN0(I-1 ,M-1 ) 

Q  -  I 

Y  -  B(I) 

IF  (P.LT.1)  GO  TO  15 
DO  10  K1-1.P 
K  -  M-K1 
Q  -  Q-1 

Y  -  Y  -  L(I,K)*X(Q) 

10  CONTINUE 
15  X(I)  -  Y*L(I,M) 

20  CONTINUE 

SOLUTION  OF  UX  -  Y 

NP1  -  N+1 
DO  AO  11-1, N 
I  -  NP1-I1 
P  -  MIN0(I1-1 ,M-1) 

Y  -  X(I) 

Q  -  I 

IF  (P.LT.1)  GO  TO  35 
DO  30  K1-1.P 
K  -  M-K1 

Q  *  Q+1 

Y  -  Y  -  L(Q,K)*X(Q) 

30  CONTINUE 
35  X(I)  -  Y*L(I,M) 

AO  CONTINUE 

RETURN 
END 
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